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Barriers to Achieving 

Textbook Multigrid Efficiency (TME) in CFD 

A chi Brandt t 


Abstract 

“Textbook multigrid efficiency” (TME) means solving a discrete PDE prob- 
lem in a computational work which is only a small (less than 10) multiple of the 
operation count in the discretized system of equations itself. As a guide to attain- 
ing this optimal performance for general CFD problems, the table below lists every 
foreseen kind of computational difficulty for achieving that goal, together with the 
possible ways for resolving that difficulty, their current state of development, and 
references. 

Included in the table are staggered and nonstaggered, conservative and non- 
conservative discretizations of viscous and inviscid, incompressible and compress- 
ible flows at various Mach numbers, as well as a simple (algebraic) turbulence 
model and comments on chemically reacting flows. The listing of associated com- 
putational barriers involves: non-alignment of streamlines or sonic characteristics 
with the grids; recirculating flows; stagnation points; discretization and relax- 
ation on and near shocks and boundaries; far-field artificial boundary conditions; 
small-scale singularities (meaning important features, such as the complete air- 
plane, which are not visible on some of the coarse grids); large grid aspect ratios; 
boundary layer resolution; and grid adaption. 

Introduction (by James L. Thomas, NASA LaRC) 

Computational fluid dynamics (CFD) is becoming a more important part of 
the complete aircraft design cycle because of the availability of faster computers 
with more memory and improved numerical algorithms. As an example, all of the 
external cruise-surface shapes of the new Boeing 777 wide-body subsonic transport 
were designed with CFD [Rl], The cruise shape of such a vehicle is designed to 
minimize viscous and shock wave losses at transonic speeds and can be analyzed 
with potential flow methods coupled with interacting boundary layers. Off-design 
performance associated with maximum lift, buffet, and flutter and the determina- 
tion of stability and control derivatives, involving unsteady separated and vortical 
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flows with stronger shock waves, are determined largely by experimental methods. 
Computational simulations of these flowfields require the use of Reynolds-averaged 
Navier-Stokes (RANS) methods; these computations for high-Reynolds flows over 
complex geometries are very expensive, the turnaround time is too long to impact 
the design cycle, and the turbulence models for separated flows have a high degree 
of variability. Thus in these areas experiments, rather than computations, are 
preferred for reasons of cost and uncertainty 

Inroads are being made into these off-design areas with RANS methods. A 
major lesson learned from industrial use of RANS methods is that both the numer- 
ics and the physics must be improved substantially for a new procedure to replace 
an older procedure. Also, there is a synergistic interplay between the speed of 
the simulation and the fidelity of the turbulence model, since a larger parameter 
variation and/or model formulation can be explored on fine enough grids with 
a faster simulation. For example, the TLNS3D Navier-Stokes code [R2] found 
its way into use because it was the first three-dimensional Navier-Stokes code to 
show true multigrid performance, in which the cost scales linearly with the num- 
ber of unknowns, and it incorporated a better turbulence model than the algebraic 
models then in use. Solutions with 1 million grid points could be converged in ap- 
proximately 1 hr of Cray-2 time, which allowed spatial convergence studies to be 
conducted to ensure that the level of truncation error is sufficiently low, and the 
prediction of the angle of attack to attain a desired lift coefficient was improved 
over interacted potential methods [R3]. The faster turnaround of the multigrid 
procedure enabled the extension and calibration of the original two-dimensional 
turbulence model to three-dimensions, thus allowing a more accurate prediction 
of the transonic shock/boundary-layer interaction. 

The current RANS solvers with multigrid require on the order of 1500 residual 
evaluations to converge the lift and drag to one percent of their final values for 
wing-body geometries near transonic cruise conditions. Complex geometry and 
complex physics simulations require many more residual evaluations to converge, 
if indeed convergence can even be attained. It is well-known for elliptic problems 
that solutions can be attained using full multigrid (FMG) processes in far fewer, on 
the order of 3-6, residual evaluations; this efficiency is known as textbook multi- 
grid efficiency (TME). Thus, there is a potential gain of two orders of magnitude 
in operation count reduction if TME could be attained for the RANS equation 
sets. This possible two order of magnitude improvement in convergence represents 
an algorithmic floor since it is unlikely that faster convergence for these nonlin- 
ear equations could be attained. This algorithmic speed-up, however, coupled 
with further increases in computational speed can open up avenues and accelerate 
progress in many areas, including: the application of steady and time-dependent 
simulations in the high-lift, off-design, and stability and control areas; the usage 
of RANS solvers in the aerodynamic and multidisciplinary design areas; and the 
development of improved turbulence models. 

The RANS equation sets are a system of coupled nonlinear equations which 



are not, even for subsonic Mach numbers, fully elliptic, but contain hyperbolic 
factors. The theory of multigrid for hyperbolic and mixed-type equations is much 
less developed than that for purely elliptic equations. Resolution of complex ge- 
ometries and the thin boundary layers at high Reynolds number cause the grid to 
be highly irregular and stretched, leading to a slowdown in convergence. Discon- 
tinuities, such as shocks and slip surfaces, introduce additional difficulties. These 
difficulties are illustrated in the sketch in Fig. 1 for a typical multi-element sec- 
tion of a three-dimensional wing with the flaps deployed at takeoff and landing 
conditions. Overcoming these difficulties poses a formidable challenge, especially 
because in order to attain optimal and robust convergence rates for the applica- 
tions of interest in aircraft design, they must all be overcome. 

Brandt, in 1984 [G84], summarized the state of the art for attaining multi- 
grid performance for fluid dynamics. Since that time, there has been considerable 
progress in the field, although optimal results have only been shown for inviscid 
flows, viscous flows at low Reynolds number, and simple geometries. The method- 
ology and theory that Brandt and others have developed is applicable to the RANS 
equations and can lead to optimal convergence rates; however, a rational and sys- 
tematic attack on the barriers which stand in the way needs to be mounted. The 
purpose of this paper is to delineate clearly the barriers which exist to attaining op- 
timal convergence rates for solutions to the fluid dynamic equations for complex 
geometries. The following sections identify the barriers, possible solutions, and 
current status of the problem. The paper is intended as a guide to attaining the 
optimal convergence goal and is written for the most part in a tabular form so that 
new solutions and updates to the current status can be made. When completed, 
the document is intended to list every type of computational difficulty encountered 
on the road to attaining TME for RANS and the solution paths taken. The in- 
sights, lessons learned, and methodologies gained from aerodynamic applications 
should be applicable to other areas such as acoustics, electromagnetics, hypersonic 
propulsion, and aerothermodynamics. 

Preliminary comments 

The table below does not refer to a vast literature on multigrid methods 
in CFD (see for example [AJ]), in which enormous improvements over previous 
(single-grid) techniques have been achieved, but without adopting the systematic 
TME approach. This ' approach insists on obtaining basically the same ideal ef- 
ficiency to every problem, by a very systematic study of each type of difficulty, 
through a carefully chosen sequence of model problems. Several fundamental tech- 
niques are typically absent in the multigrid codes that have not adopted the TME 
strategy. Most important, those codes fail to decompose the solution process into 
separate treatments of each factor of the PDE principal determinant, and there- 
fore do not identify, let alone treat, the separate obstacles associated with each 
such factor. Indeed, depending on flow conditions, each of those factors may have 
different ellipticity measures (some are uniformly elliptic, others are non elliptic 
at some or all of the relevant scales) and/or different set of characteristic surfaces, 
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requiring different combinations of relaxation and coarsening procedures. 

The table deals only with steady-state flows and their direct multigrid solvers, 
i.e., not through pseudo-time marching. Time-accurate solvers for genuine time- 
dependent flow problems are in principle simpler to develop than their steady- 
state counterparts. Using semi implicit or fully implicit discretizations, large and 
adaptable time steps can be used, and parallel processing across space and time is 
feasible [R88]. The resulting system of equations (i.e., the system to be solved at 
each time step) is much easier than the steady-state system because it has better 
ellipticity measures (due to the time term), it does not involve the difficulties 
associated with recirculations, and it comes with a good first approximation (from 
the previous time step). A simple multigrid “F cycle” at each time step can 
solve the equations much below the discretization errors of that step [Par]. It is 
thus believed that fully efficient multigrid methods for the steady-state equations 
will also yield fully efficient and highly parallelizable methods for time-accurate 
integrations. 
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Textbook Multigrid Barriers 
(After Brandt, 1 997) 


Far-Field Boundary Conditions Grid Aligned with Characteristics 


Grid Skewness 


Shock 


Inviscid Transonic 
Flows and Sonic 
Lines 


High-Aspect-Ratio Grids 


Boundary-Layer Separations 


Grid Near Boundaries 


Laminar Separation 
bubbles 



Boundary-Layer 

Separations 



Attaining Ideal Multigrid Efficiency in CFD: Difficulties and Cures 


Difficulty 

(•) Uniformly elliptic scalar 
equation on uniform grids 
in general domains 

(•) Nonlinearity 

(•) Fluid dynamics - general 


(•) Non-scalar PDE systems 


Possible Solutions 

Multigrid cycles, guided by local mode analyses 

+ FMG 
FAS 

+ FMG continuation 

See a review in [R88, §2]; at some points it is not 

fully up to date, but it concisely summarizes the 

main procedures needed for obtaining TME 

(1) General rules for the inter-grid transfers are 
given in [G, §4.3], with some more details in 
[RQMA, §3.3] 

(2) General approach to the design of relaxation, 
based on the operator principal matrix L and 
on the factors of det L. Distribution matrix 
M and weighting (or preconditioning) matrix 
P are constructed so that PLM is triangular, 
containing the factors of det L on the main 
diagonal (separated from each other as much 
as possible, to avoid the complication descri- 
bed next). This (if necessary - together with 
the technique described next), leads to decom- 
posing relaxation into simple schemes for the 
(scalar) factors of det L 

(3) For systems of PDE which are of mixed type 
(elliptic-hyperbolic) another possibility is to 
introduce new unknowns in terms of which 
elliptic and hyperbolic parts are separated 


Status 

TME demonstrated 1971 [B73], 
[B77] and rigorously proved 
[RLMA], [RQMA] 

Demonstrated 1975 [South], [B77], 
[G, §8.3.2] 


TME demonstrated in a number 
of cases (see below). TME proved 
for uniformly elliptic systems 
[RLMA], [RQMA] 


TME demonstrated for incom- 
pressible and compresisble 
cases [T1]-[T5] 


Difficulty 


Possible Solutions 


Status 


• Product operator: an equation 
LU = /, where L = L2L1, with 

the Fourier symbols 
Lj( 0 ) = e -w-x/h Le ie-x/h 

Assume a relaxation process for 
Lj is given, with the amplification 
factor Pj(0) and the smoothing 
factor Jlj, ( j = 1,2) 


Two possible approaches: 

(1) Introduce an explicit new unknown function V, 
replacing the equation with the pair of equations 
L\U — V = 0 and L 2 V — /, throughout the MG 
solution process (including, e.g., transferring 
residuals of both equations to coarse grids and 
correcting both u and v by interpolations from 
the corresponding coarse-grid values). The smoo- 
thing factor for this process is Ji = maxf/U, /7 2 ) 

(2) Use V only as an auxiliary function in relaxation. 
That is: starting with v = Zqu, where u is the 
current approximation to JJ , perform V 2 sweeps 
on the equation L 2 V = /, yielding a new value v. 
Then perform v\ sweeps on the equation L\u = v. 
The resulting amplification factor is 

m = w(*r + [1 - w(«r]£i(0)- w, 

so in scalar cases p < p'y + /j, 2 ^ 2 


Not tried 


Difficulty 


Possible Solutions 


(•) Smoothing for special CFD 
systems 

• Cauchy Riemann on staggered grid 

d x d v 


L = 


dy -d x 


Stokes on staggered grid 


L = 


• Stokes, non-staggered 

(1) Quasi-elliptic discretization 

/-A 0 


L = 


0 -A 

V d 2h dy h 


df 


0 


with averaging of the resulting 
pressure 

(2) /i-elliptic discretization, e.g. 


L = 0 


-A 

0 

9f \ 

0 

-A 

df 

dl h 

q2Ii 

Uy 

—ah 2 A / 


M = distribution operator 
P = preconditioner 

(1) M = L, P = I 

(2) P = L, M = I 


/I 0 -d x 


-du 


/-A 

0 

d x \ 

(1) M = 

0 1 

0 

-A 

dy 


\0 0 

d x 

Dy 

0 / 

/ 

/ 1 

0 




(2) P = 

0 

1 


P = I 


0 

0 


M = 1 


\d x d, 


y 


Analogous to the staggered case: 

1 0 ~dl h 
M = ( 0 1 

.0 0 -A 


-df 


No modifications of the FMG 
algorithm is required, even in the 
quasi-elliptic case (as explained in 
[G84, §18.6]). In generalization to 
NS, pressure averaging is required 
of coarse-level results before their 
interpolation to the next finer level 
(whenever the coarse-level employs 
the quasi-elliptic discretization) 


Status 

(1) TME demonstrated [BD], [Dinar] 

(2) TME validated [T6] 

(1) TME demonstrated [BD], [Dinar] 

(2) TME validated [T6] 

(1) In a quasi-elliptic approach, TME 
demonstrated [G84, §18.6], [quasi] 

(2) TME demonstrated? (perhaps by 
J. Linden) 



Difficulty 


Possible Solutions 


Status 


• Non- conservative incompre- 

ssible Euler 

/ U. ■ V 0 d x \ 

L = I 0 u • V d y 

V d x d y 0 / 

(similarly 3D) on staggered 
grid, second (or higher) order 
discretization 


• Low- Reynolds Incompre- 
sibble NS, staggered or not 

• High- Reynolds Incompre- 
ssible NS, staggered or not 


( 1 ) Employ cycle index 7 = 2 P , where p 
is the order of discretization 



/I 

0 

dx 

(2) M = 

0 

1 

-dy 


Vo 

0 

u ■ V 


For each of the momentum equations 
employ a relaxation scheme which is 
fast converging for the advection 
operator u • V (i.e., converging fast 
not only for h-f , but also for smooth 
characteristic components; see discu- 
ssion of advection below) 

(3) Use canonical variable (u,v,P) on sta- 
ggered grid, where P — (u 2 + y 2 )/2 + p. 
Upwind only P, use central discretiza- 
tion for (u,v). Relaxation is marching 
for P , and weighted (preconditioning) 
for (u, v ) 

Fully analogous to Stokes solvers: just 
replace A in L by Q — — R -1 A + m • V 

Fully analogous to Incompressible Euler 
(outside boundary layers: see discussion 
on such layers below): just replace u ■ V 
everywhere with Q 


(1) TME for first-order discretization 
using W cycles shown in [BD], [Dinar] 

(2) TME demonstrated for 2D entering 
flows with second-order discretization 
[BY2] and for recirculating flows with 
first-order discretization [BY3] 


(3) TME in [T1-T3] 


TME demonstrated 1978 [BD], [Dinar] 

TME demonstrated for first-order discreti- 
zation on staggered ([BD], [Dinar] and 
non-staggered grids [G84, §19.5], and for 
second-order staggered discretization [BY2] 



Difficulty 


Compressible Euler , non-conservative, on staggered grid: 
The subprincipal operator on (uj, 112 , U 3 , p, £,p) is 



/ pu. • v 

0 

0 

0 

0 

di\ 


0 

pu ■ V 

0 

0 

0 

d 2 

T _ 

0 

0 

pu • V 

0 

0 

d 3 

Jj — 

P 2 d 1 

P 2 d2 

P 2 d^ 

pu ■ V 

0 

0 


pd 1 

P&2 

pd 3 

0 

pu ■ V 

0 


V 0 

0 

0 

—dp/ dp 

-dp/de 

1 / 

det L 

= p 5 (m. • 

Y) 3 ((n. 

• V) 2 - 

a 2 A) 



2 

a = 

dp P_ 

dp p 2 

dp 

Ik' 

M o = | 

\u\/ a 




p, e, p defined at cell centers, 

U j - at center of cell faces perpendicular to the i-th coor- 
dinate 


2D Compressible Euler, nonconservative and conserva- 
tive, staggered grid, using canonical variables (u,v, S , H ). 
Structured and unstructured grids 


2D/3D incompressible and compressible Euler: Canonical 
variables in which velocities are replaced by vector poten- 
tial representation. Nonstaggered structured and unstruc- 
tured grid 


Possible Solutions 


M = 


Z 1 

0 

0 

0 

0 

Vo 


0 0 
1 0 
0 1 
0 0 
0 0 
0 0 


0 0 
0 0 
0 0 
1 0 
0 1 
0 0 


-p(u-V)di \ 
-p(u ■ Y)^2 
-p(u • V)d 3 
-p 2 A 
—p A 

p 2 {u • V) 2 / 


The advection and full-potential operators 
are each relaxed by one of the approaches 
described for them below. (The semi 
coarsening described there would then be 
used as an inner multigrid cycle for 
relaxing one factor of the determinant, to 
be distinguished from the outer multigrid 
cycle, which can use full coarsening.) 


Use (u,v) at cell edges, H at middle of cell, 
S at vertices. Upwind only S at momentum 
equations. Relax S, H by marching. (u,v) 
by a. weighting relaxation. Crocco’s form is 
used here to define relaxation 

All variables at cell nodes. Relax hyperbo- 
lic quantities using marching. Relax vector 
potential using point Gauss-Seidel 


Status 

Not tried 


TME in [T2-T5] 


TME acheived 
(unpublished) 
for interior 
and exterior 
flows in 2D, 
interior in 3D 



Difficulty 


• Compressible Navier- Stokes, non-conservative. 

The subprincipal operator on (u\, w 2 , «3, p, £,p) is 

Qn — Adn —Xdi2 — A<9i3 0 0 

-Xd -21 Q/j — Ad 22 — A ^23 0 0 

-A 031 -A ^32 Qh - A ^33 0 0 

p 2 d 1 P 2 d '2 P 2 d 3 Qo 0 

pd l pd 2 pd 3 0 Q k 

0 0 0 —dp/ dp —dp/de 

where Q a = -aA + pu • V, A = A + p, A = §/i, 

K — kjc v (coefficient of thermal conductivity divided 
by the specific heat at constant volume), 

■"J det L s = Qj t det L c , where L c is the “core operator” 


/ Qo 

0 

—p 2 A 

L c =l 0 

Q k 

—p A 

\ -dp/ dp 

- -dp/de 

Q n+\ 


At standard conditions of laminar air flow the 
Prandtl number 77 /,/ k 0 . 72 ; for turbulence 7 /.// k ps 0 . 9 , 
with 7 = Cp/c v = 1.4 



Possible Solutions Status 

(1) Where A, p, k < ph\u\ relax as in Euler above Not tried 

(2) Otherwise use 

1 0 0 0 0 -di 

0 1 0 0 0 -d 2 

0 0 1 0 0 -d 3 

0 0 0 1 0 0 

0 0 0 0 1 0 

Adi A d 2 Ad 3 0 0 Q^ + j 

relaxing each Q fl by one of the approaches 
described for the ad vection- diffusion below, 
and L c by procedures discussed for it below 
(in the chapter on non-elliptic operators) 





Difficulty 


Status 


Non-conservative not staggered 
Euler and NS 


Conservative discretization of 
any of the above systems 


Possible Solutions 

(1) Probably similar to the staggered (cf. 
transition from staggered to non-stag- 
gered in Stokes) 

(2) In the 2D incompressible case: 
Premultiply L by a projection operator 
P, obtaining a Poisson equation for the 
pressure. Solve pressure equation with 
multigrid and the advection equation by 
marching downstream. 

Apply a prefactor P such that PL has 
principally the above non-conservative 
form. See, however, the difficulty 
associated with FDA decompos ability 
(discussed in the chapter on non- 
elliptic operators), which may arise 
with such PL operators 


TME demonstrated for 2D 
incompressible Euler [RSS] 
in the cases of channel 
(with bump) and airfoil flows 


Mentioned in [G, §3.4], but 
not tested 



Difficulty 


Possible Solutions 


Status 


® Non-elliptic operators, or more pre- The DGS relaxation of the full flow 

cisely: small ellipticity measures at some equations allows a specific individual 
(e.g., large) scales. The main operators treatment for each of these cases, 

of interest here are taking into account its particular set 

(1) The advection operator (or, similarly, of characteristic 
the convection-diffusion operator at 

large Reynolds numbers). 

(2) The near-sonic full-potential opera- 
tor or more generally the core opera- 
tor L c . 

(See below a separate discussion of aniso- 
tropies caused by the discretization) 

• Grid aligned with the characteristics Block (e.g., line or plane) or ILU TME demonstrated in 

relaxation schemes and/or semi- many cases 

coarsening, possibly in alternating 
directions, guided by mode analyses 
[B77], [Stages] 

• Distinguishing different regimes Running separately the relaxation subroutine of a 

(open vs. closed characteristics) given non-elliptic factor can 

(1) Separately check its convergence properties 

(2) Produce a scalar a fa 1 at regions of open charac- 
teristics and cr <C 1 on closed characteristics 

(such as separated flow zones) 



Difficulty 


Possible Solutions 


Status 


• Non-aligned grids, with open 
characteristics (e.g., entering flow): 
The main difficulty is the shorter 
distance (along the characteristics) 
for which a coarser grid still appro- 
ximates some smooth solution compo- 
nents (characteristic components 
with intermediate cross-characteris- 
tic smoothness) [NESP], [BY1] 


Three possible approaches, all guided by half-space 
two-level FMG mode analysis, using for simplicity 
the first Differential Approximation (FDA) to the 
discrete operator [NESP], [G, §7.5]: 

(1) Downstream-orderecl relaxation marching 
[R88, §2.3]. (Suitable only for the advection 
factor, sometimes still requires W cycles, and 
not very good for massively parallel 
processing). In the case of an 0(h p ) discretization 
which is not purely upstreamed, relaxation 
should involve a predictor-corrector downstream 
marching. If the predictor order is q , the corrector 
should be applied at least p/q times. 

(2) Semi-coarsening (especially for the near-sonic 
full-potential operator). Better suits massive 
parallel processing 

(3) Cycle index = 2 p / m , where p is the order of 
discretization and m is the order of the differential 
factor. (Suitable actually only for the advection 
operator, for which m — 1; especially attractive 

p — 2 in 3D; not requiring ordered relaxation, but 
still disadvantageous for massively parallel proce- 
ssing because of the high cycle index) 


(1) TME demonstrated in 
[BY2], and in Ruge’s 
recent calculations, 
both for the advection 
operator by itself and as 
part of the incompre- 
ssible Euler system 

(2) TME has been shown for 
the sonic full-potential 
operator [sonic] 

(3) For p = 1, TME has - 
been shown on various 
occasions. For p = 2, 
should be tried 



Difficulty 


Possible Solutions 


Status 


• The mixed convection- diffusion opera- 
tor with order p approximation, having 
natural viscosity v and artificial visco- 
sity ah p 

• Closed characteristics (recirculating 
flows). Here uniformity of viscosity 
(including numerical viscosity) is 
important for accuracy. The size 

of viscosity is less important here 
(except at resolved boundary layers, 
discussed below). In fact, a uniform 
0(h) artificial viscosity can yield 
higher order approximations. Full 
convergence may also be less impor- 
tant here (steady state may take 
exceedingly long to be attained in 
reality, if at all) 


• Full-potential operator ( u ■ V) 2 — a 2 A, 
M 0 = |«|/a ^ .7 (uniformly elliptic) 


Treatment as elliptic operator on levels where 
v > (2 P • 4 — 5 )ahP and as the non-elliptic 
advection operator otherwise 

Using the above-mentioned scalar a, form a 
<7-dependent convergence test, to tell between 
slowness of open and closed characteristics 
(and possibly ignore the latter). Also based 
on <7, at recirculation regions use uniform 
(explicit) 0(h) numerical viscosity, with con- 
tinuation from large to small viscosity integra- 
ted into the FMG algorithm. The cycles can 
employ one of the following 3 options. 

(1) DCW method (using Defect Corrections 
within W cycles), with suitable over- 
weighting of residuals [BY3]. Suitable only 
for 0(h) discretizations. 

(2) Effectively downstream relaxation ordering 
(using alternate- direction sweeps) and 
doubling of transferred residuals (for 0(h) 
discretization) [YVB], 

(3) Semicoarsening, generally similar to [sonic] 

Any classical algorithm is suitable, but the 
algorithm of the next case is also adequate 


Not precisely tried 


TME cycles by methods 
(1) and (2) were shown 
in [BY3] and [YVB] 
respectively. Method (3), 
which should be best for 
massive parallelization, 
has not been implemented 


TME well established 



Difficulty 


Possible Solutions 


Status 


Full potential .7 £ Mq < 1.4 


Relaxation marching downstream (for 
transition to the supersonic case below) 
together with semicoarsening in the 
characteristic (cross-stream) direction 


TME shown for the 
case M = 1 [Sonic], 
Other cases have not 
yet been implemented 


Full potential 1.4 < Mq 
(uniformly hyperbolic, with the 
stream as the time-like direction, 
and with 0(1) “Courant number”.) 


Marching in the stream direction, possibly Not yet tried? 
with a predictor-corrector procedure. 

For full massive parallelization, however, 
wave methods (extending standing wave 
methods [Ira]) should be used 



Difficulty 


Possible Solutions 


• The “core operator” 



Qo 

0 -p 2 A' 

L c = 

0 

Q k —pA 


_ —dp/ dp 

-dp/de 


should be relaxed as part of relaxing 
the compressible NS system, in the 
case that p\u\h < max(A, p, «). 

In the case of alignment between the 
grid and the flow, with meshsize h\ 
and h 2 in the stream and cross- 
stream directions, respectively, and 
/i-2 5: hi (e.g., in boundary layers), 
the case where L c need be relaxed is 
when p[u\h2 < h\ max(A, p, k). 

In aerodynamics, A, p and k are 
comparable, so the case of interest is 
\u\h 2 < vh\, where v = p/ p 


Best relaxation scheme depends on the flow parameters. 

For example: 

(1) If k <C p\u\h, then Q K & Qq (in principal terms) and one 
can use DGS with 

( 1 0 p 2 A\ 

0 1 pA I 

0 0 Q 0 / 

resulting in the need to relax the first two equations 
each on an advection operator (see methods above), and 
the third equation on the operator QqQ — p 2 a 2 A. 

In the case of interest the principal part of the latter 
is [(/i + \)Qo + p 2 a 2 ] A, so it can be relaxed 
by the general method for relaxing a product opertor 
(see L = L-jLi above). 

(2) In the aerodynamics and aligned case of interest, the term 
Q |_y in L c is not principal. Therefore relaxation can easily 
be conducted with the weighting (preconditioning) matrix 

/I 0 0 

P = j — p p 2 0 

\ 0 0 1 

and the distribution matrix 

(0 1 0 \ 

M = \ 0 -pp/pe 1 

\1 0 0 / 

yielding PLM whose principal part is its main diagonal, on 
which separately appear the Laplace operator A, the 
convection-diffusion operator Qk where k = p p p 2 /(2pp e ) 

= 1.25k (for air), and a free function 


Status 


Not tried 


Not tried 



Difficulty 


Possible Solutions 


Status 


• FDA factorizability question: The decom- 

position of a system relaxation into its 
scalar factors depends on the equality of 
the different occurrences of the advection- 
diffusion operator Q (or appearing 

in PL, the prefactoring by P of a conser- 
vative discretization L. However, in rela- 
xing a non-elliptic discrete operator, impor- 
tant is not only the differential operator it 
approximates, but also its First Differential 
Approximation (FDA) terms in non-charac- 
teristic directions; e.g., the cross-stream 
numerical viscosity of Q. This may not be 
the same in the different occurenes of Q, 
putting the factorization into question 

• High order discretization (away from 
shocks) 


(1) Examining several examples of con- 
servative discretization of transsonic 
flows, the FDA terms in various 
occurrences of Q fl+ j turn out suffi- 
ciently close to each other (e.g., only 

(4% discrepancy) to allow full effi- 
ciency of the proposed relaxation 
schemes. 

(2) Conservative schemes may be desig- 
ned so that the various FDAs of 

are identical, or at least so 
that the scheme is still factori- 
zable. 


(1) “Double discretization” schemes: 

Use high-order only in calculating 
residuals transferred to the coarse 
grid, not in relaxation (unless the 
high order scheme is preferable also 
for h-f modes). 

(2) However, in relaxing non-elliptic 
factors (e.g., downstream relaxation 
marching for convection operator) 
the high order must be used (e.g., by 
a predictor-corrector downstream 
relaxation) 


Further examination is needed 


Some “genuinely multidimen- 
sional upwind” schemes turn 
out to yield factorizable 
schemes, e.g., in the subsonic 
case in the control- volume 
structured-grid context [DS2]. 
Further studies are in progress 

Introduced 1978 [BD]. Success- 
fully implemented in various 
elliptic cases (see description 
and refs in [G, §10.2]). Methods 
for non-elliptic have not been 
tested beyond second order. 

Comment: High order approxi- 
mations on unstructured grids 
are very expensive 



(1) A set of N continuity equations, 


volume is large compared with max(/? ~D;. h 1 pi\u\), 


for a simple model 


Difficulty 

(•) Shocks 

• Shock displacements associated 
with corrections from a coarse grid 
that does not resolve the shock 

• Poor h-ellipticity of high-resolution 
schemes 

• Relaxation near strong shocks 


Possible Solutions 


Obtained by a. conservative fine- to- coarse 
residual transfer plus local post-relaxation 
passes near the shock 

Construction of new, genuinely multidimen- 
sional upwind schemes 


Switching to general robust schemes (e.g., 
box Kacmarz), adding extra local passes (c.f. 
relaxation near boundaries) 


Status 


Full efficiency shown [DS] 


Developed in the context 
of unstructured triangular 
grids [DSl] 

Not tried? 



UliVy ^ L-OOl U1C KJL lliLAJllipi 

N a.vi Pt'- *-i fokf^Q pniiftf.inTiG arlrli-no* f<~* fl-ip* 


Difficulty 

Relaxation at and near Boundaries : 
Difficulties: 

(1) There is no smoothing analysis in 
case the boundaries are not aligned 
with the grid. 

(2) The fine-to-coarse residual weigh- 
ting near boundaries is generally 
very imprecise, hence the residuals 
should be reduced there more than 
in the interior. 

(3) Larger residuals are created near 
polations. 

Boundary layers (in the case that they 
need be resolved). (See also grid 
adaptation below.) 


ilDIC J 
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be demonstrated 

t — _ i - 


Possible Solutions 

A general-type robust relaxation scheme, 
e.g., box Kacmarz , throughout several- 
meshsize-wide zone near the boundary. 
Boxes size in each direction should be 
several meshsizes and the boxes should 
have substantial overlap. One can afford 
several passes of such a relaxation per 
each full interior sweep since the zone 
width is 0(h }~ £ ), with 0 < £ < 1. In 
particular, add near-boundary relaxation 
passes after the FMG interpolation 
(allowing the latter to be of lower order 
near the boundary) 

Resolved by boundary-fitted local grid 
patches, with local semi refinements: finer 
levels, in narrower layers near the boun- 
dary, have smaller cross-layer meshsizes, 
allowing the physical cross-stream visco- 
sity to dominate over the numerical one. 
Additional terms in the governing equa- 
tions (NS instead of Euler, or turbulent 
modelling, etc.) may be used in these 
patches. Downstream marching relaxation 
and cross-stream semi coarsening in the 
multigrid cycles, employed in a “A-FMG” 
kind of algorithm [G, §9.6], so that coarse 
FMG stages already include local semi- 
refinements at the boundary, thus effectively 
incorporating continuation in R e into the 
FMG stages 


Status 

For uniformly elliptic equations it 
has been proved [RLMA], [RQMA] 
and demonstrated computationally 
(for cases of reentrant corners [Bai]) 
that the interior efficiency as predic- 
ted by mode analysis (implying TME) 
can always be obtained. TME demon- 
strated (by Ruge & Brandt) for incom- 
pressible Euler on staggered cartesian 
cartesian grids 


Description in [R88, §2.4]; not imple- 
mented. The local refinement techni- 
ques for Poisson equation, with TME, 
are demonstrated in [Bai] 


Difficulty 


Possible Solutions 


Status 


Far-field artificial boundary 
conditions: requiring in some 
cases non-local absorbing 
boundary conditions (ABC) for 
some wave factor. [Has any 
MPer had experience with this 
difficulty?] 


Small-scale singularities invisible 
on the next coarser grid, such as 
small “islands” or “holes” in the 
domain (e.g., an airplane smaller 
than the meshsize of some coarser 
grid) or small BC features (e.g., 
small regions of Neumann BC in 
otherwise Dirichlet BC) 


Increasingly coarser grids covering increasingly 
larger domains. The size of each domain is 
based on accuracy-to-work. Optimization 
criteria (similar to those in [B77, §8], [G, §9.5], 
implying also a natural criterion for the largest 
needed domain. On interior boundaries (boun- 
daries of a grid residing in the interior of the 
next coarser grid) the solution is interpolated 
from the coarser grid. On such boundaries, if 
ABC is at all needed, only high-frequency 
components need be absorbed, for which the 
ABC are local , and can be enforced as part of 
the relaxation process (of the corresponding 
wave factor) 

Local relaxation passes around the singularities 
after return from the next coarser grid, together 
with either one of the following three devices: 

(a) Enlarging the singularity on the coarser 
grid. 

(b) Modifying the interior coarse-grid equation 
near the singularity. 

(c) If the coarse grid equations are not modi- 
fied, convergence is slow, but slow to con- 
verge are just few very special components. 
Hence slowness can be eliminated by recom- 
bining iterants. 


Details of the algorithm have been 
worked out, and TME (or its equi- 
valent accuracy-to-work relation) 
has recently been demonstrated (by 
Brandt & Danowitz) for the 2D 
Poisson equation in the unbounded 
plane. Techniques for non-elliptic 
or indefinite cases have not been 
systematically studied. 


TME shown in elliptic cases [Rec] 



Difficulty 


Possible Solutions 


Status 


(•) Grid-induced slow 


• Large aspect ratios 


• Expanding grids 


convergence One can avoid many of the following maladies 
by using suitable multi-grid structures (descri- 
bed below under “grid adaptation”) 

Either one of the following: 

(1) Block (part-line or part-plane) relaxation, 
analyzed by mode analysis [B77]. 

(2) Semi coarsening [Arl], [Stages, §3.2] (often 
natural, since the large aspect ratio is in 
the first place created by semi refinements) 
with relaxation “semi smoothing” analysis 
[Stages, §2.1], [G, §3.3], 

(3) Combinations of block relaxation in some 
directions and semi coarsening in others 

Relaxation marching in the direction of increa- 
sing meshsize; or distributive relaxation [Njm, 
§ 6 ] 


TME has been shown in 
a variety of elliptic cases 



Difficulty 


Possible Solutions 


Status 


• Grid adaptation 


(•) Stagnation point (causing an 
instability in the coarse-grid 
corrections) 


Use local multigrid levels in creating any 
desired local refinement, aspect ratio, 
boundary fitting or even flow fitting (see 
[R88, §2.7]). Base refinement criteria on 
the fine-to- coarse multigrid correction (r). 
Adaptation can be integrated into the 
A-FMG algorithm together with proper 
(e.g., Reynolds-number) continuations 

Coarse-grid numerical viscosity depending 
on the average (e.g., “full- weighting”) of 
the fine-grid numerical viscosity (not on 
its injected value) [BY3, §4.5] 


Introduced in [B77] and [G], but 
tried only for Poisson equation 
near singularities [Bai] 


TME shown in an example [BY3] 
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